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Abstract 

We study the effects of generalised surface disorder on the mon- 
omer-monomer model of heterogeneous catalysis, where disorder is 
implemented by allowing different adsorption rates for each lattice 
site. By mapping the system in the reaction-controlled limit onto a 
kinetic Ising model, we derive the rate equations for the one and two- 
spin correlation functions. There is good agreement between these 
equations and numerical simulations. We then study the inclusion 
of desorption of monomers from the substrate, first by both species 
and then by just one, and find exact time-dependent solutions for the 
one-spin correlation functions. 
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I. INTRODUCTION 

Diffusionless surface-reaction models were first introduced by Ziff, Gulari 
and Barshad who investigated a monomer-dimer reaction corresponding 
to the chemical reaction 2C0 + 02^ 2CO2 on a catalytic surface. A well- 
studied variant 0, 0, Q employs the simpler monomer-monomer reaction, 
described by 



A -I- ^ A 

-'^gas ^ ' ■'^surface 

f) _|_ Q r) 

-'-'gas ~r O J^surface 

^surface ~l~ -^swr/ace ^ ^-^gas ~l~ 2S', (l) 

where S denotes an empty site. This process exhibits a kinetic phase when 
there are equal propensities of A and B species, in which the long-time ki- 
netics become dominated by domain coarsening. Mean- field analysis f^, in 
which every site is taken to be connected to every other site in a 'complete 
graph', demonstrated that finite lattices will always saturate - that is, the 
lattice will either become full of A's, or full of B's, and the process will stop. 
Krapivsky ^ recently solved the model exactly in the reaction-controlled 
limit /c/j — i> cx) by mapping the system onto the standard Ising model. 

Many enhancements to these models have been studied with a view to 
more closely modelling actual chemical processes, including nearest neigh- 
bour excluded adsorption [0 and surface diffusion |^. However, only re- 



cently have the effects of surface disorder been touched upon by Frachebourg 
et.al. 0. They chose to model a disordered surface by taking a lattice of 
two different types of site, one which favours adsorption by the A-species, 
and one which favours adsorption by the B's. They showed numerically that 
such disorder allows for a reactive equilibrium in two dimensions. 

In this paper, we extend the analytical method used in |^ to a general 
form of surface disorder, based on but allowing for a range of different 
types of site in the lattice. Furthermore, we also investigate separately the 
effects of desorption in the system. All the results presented are for the 
physically relevant case of two dimensions. 

This paper is organised as follows. In section II we define the model 
and derive the general rate equations for the n-spin correlation functions. 
In section III, these equations are applied to a model similar to that in 
and their solutions are compared to numerical simulations. In sections IV 
and V we include the effects of desorption, first by both species and then 
by just one, and derive exact solutions. The conclusions are summarised in 
section VI. 

II. RATE EQUATIONS 

We consider the surface reaction A + B ^ 2S on a periodic L x L square 
lattice, ignoring the effects of diffusion and desorption. For simplicity, we 
take the reaction-controlled limit, where the adsorption of A and B species is 
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taken to be infinitely fast so tliat tlie substrate is always full. The algorithm 
employed here is to select a nearest-neighbour (NN) pair at random, check 
for an AB-reaction, and, if so, remove the particles and immediately refill 
both sites. 

With the usual homogeneous model, the probability of filling a site with 
and A or B is independant of the site chosen - in this model, however, that 
probability is allowed to vary. Specifically, we introduce the site inhomo- 
geneity matrix Pij, < Pij < 1 Vz, j, such that the probability of filling the 
site with an A is given by Pij (or, equivalently, a probability 1 — P^j of 
filling the site with a B). 

Since in the reaction-controlled limit each site has only two possible 
states, we can map this model onto an Ising model with mixed Glauber- 
Kawasaki dynamics identifying A's with Sij = +1 and B's with Sij = —1. 
The master equation for P{S,t), the probability distribution for the system 
to be in the state S = {Sij} at time t, is 

fp{S,t) = Y.[U^AF^JS)P{F,,S,t)-U.,{S)P{S,t)] 

+ E[ ^AF^,F^-^l,S)P{F,,F,^^,S, t) - V,,{S)P{S, t) ] 

id 

+ J2[WniF^^F,^+iS)PiF,^F,,^,S,t)-W,,iS)PiS,t)]. (2) 

id 

The flip operator Fij acts on the system state vector S by flipping the 
sign of the Sij component, leaving the remaining components unchanged. 
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Uij corresponds to Glauber spin-flip dynamics |T^, whereas Vij and Wij cor- 
respond to Kawasaki exchange dynamics. Equation (H) is identical the the 
homogeneous case, except that now the full expressions for Uij, Vij and Wij 
are given by 



4:TiUij = {1 - SijSi+ij){l - dij + Sij{l - afj)} 

+ (1 - SijSi^ij){l - di^ij + ^^^(l - af_^j)} 

+ (1 - SijSij - Cij + Sij{l - 6+)} 

+ (1 - - e,,_i + S,,il - b±_,)}, (3) 

'^T2Vij = {1 - SijSi+ij) [dij + ar.Sij] , (4) 
4:T2Wij = {1 - SijSij+i) {^Cij + bijSij^ , (5) 

where the constant coefficients dij and Cij are related to the inhomo- 

geneity matrix Pij, 



afj — Pi+ij ± Pij, 



rl — P J- P — op P 

Uiij 1 ij -r 1 i+lj jj'-' 



Cij — Pij + Pij+l — 2PijPij^i. (6) 

We proceed by deriving the rate equations for the one and two-spin cor- 
relation functions, where the general n-spin function is given by 



i^hjl ■ ■ ■ ^injn) — X! ^hjl ■ ■ ■ Sinj„P{S, t) . (7) 
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Using this and some lengthy but straightforward calculations result 
in the following hierarchy of differential equations, using the renormalised 
time scale r defined by = rf"*^ + T2^, and setting ri = T2, 



^^^(•^^j) = \J{S^J) + il~2P,J){S,,{A,,S,,}), (8) 
4:T—{SijSki) = {Aij + Aki){SijSki) 

+il-2Pi,){S,,SM{A,^S,,}) 
+{l-2Pki){S,,SM{AMSM}). (9) 
... for \i- k\ + \j - /| > 1 

Here, Aij{Sij) = -4:{Sij) + {Si+ij) + {Si-ij) + {Sij+i) + {Sij^i) is the 
discrete Laplacian. For |2 — A;| + \ j — /| = 1, i.e. for nearest neighbour 2-point 
correlations, the rate equation has a more complex form. For example. 



6 



~\~{Si-ijSi^ij) + {SijSi^2j) + {SijSi^ij^i) 

+ (S'jjS'j+ij_i) + {Sij-iSi^ij) + {Sij^iSi-^-ij) 

3 



+(1 — 2Pij) i^{Si-ijSijSi+ij) — -(Sj+ij) 

+ (1 - 2Pj+ij) |(5'ijS'i+ijS'i+2j) ~ 2'^*^*-''^} 

+2(1 -2di,). (10) 

In the homogeneous hmit Pij ^, the results in are recovered. 

III. TWO-SITE DISORDER 

We now turn to the case where Py can take just two different values, p 
or q = 1 — p, with an equal number of p-sites and g-sites. This corresponds 
to the model given in ^ with equal fluxes of A and B species, e = |p — || 
and c_ = c+ = |, using the notation given there. 

Since p + q = 1, the global dynamics of the system must be unchanged 
under the transformation {p,q) (1 — p, 1 — g) = {q,p)- This symmetry 
means that the system cannot favour one state over the other, and so the 
average of (5.^) taken over the entire L x L lattice, j2 Y^ij {Si j) , will always 
tend to zero in the L — 00 limit. An important consequence of this is that 
if a finite system always saturates, then it does so with equal probability 

7 



of saturating either to every site being +1, or every site being -1, and so 
{Sij)\t=oo = Vi, j, regardless of whatever Pij may be. If a reactive steady- 
state occurs - that is, if the average saturation time diverges at least as fast 
as Pi - then it should be expected that (Sij) may be non-zero for t — > oo 
(if p 7^ |). It is the purpose of this section to apply the rate equations derived 
in section II to predict the equilibrium value of {Sij) on p-sites in any such 
non-trivial steady state. 

Although the concentrations of p-sites and g-sites are equal, different 
arrangements of the sites can dramatically alter the long-time dynamics of 
the system. For instance, choosing to split the lattice into two alternating 
c(2 X 2) sublattices, with one sublattice full of p-sites and the other full of 
g-sites, results in a system with no non-trivial steady states for p 7^ or 1. 
Since saturation always occurs, {Sij)t=oo = on either type of site. 

A more informative model can be constructed by randomly arranging the 
p and q sites. This allows for regions of p-sites, which will all tend to be fixed 
into the same state, and regions of g-sites, which will all tend to be fixed into 
the other state, to 'pin' the dynamics into a reactive equilibrium. Although 
exact analysis of this model is obviously impossible, a useful approximation 
can be made by assuming that every site is surrounded by exactly 2 p-sites 
and 2 g-sites. It is then possible to write down (H) and ([l^) for the two sorts 
of site, {Sij)p and {Sij)q, and the various two-point functions. 

To obtain a closed set of equations, further approximations must be made 
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to reduce the three-point functions in ( |T0| ) to one and two point functions. 
The obvious choice is 



{SijSklSmn) ~ {SijSkl){SklSmn), (H) 

but this has the unwanted side-effect that {Sij)p + {Sij)q 7^ 0, something 
which cannot be true since p + q = 1. To restore the required symmetry we 
must also include an alternative three-point approximation, 

{SijSklSmn) ~ {Sij) {SfclSmn) ■ (l^) 

For greater clarity, we denote the one-spin correlation function 
{Sij)p = —{Sij)q by Hp, the two-spin correlation function between two NN 
p-sites (or, equivalently, two NN g-sites) by Zpp, and use Zpg for the two-point 
function between nearest neighbour p and q sites. Setting r = 1, we can now 
obtain a closed set of equations. 



2^yp = -2yp+ {l-2p){zpg + Zpp-2), (13) 

{(1 - 2p)yp + 3zpp}{zpp + Zpg) + (2 - Apq), (14) 
4^2;pg = 4:pq- {Apq + 6)zpg + 3{l-2p)yp + 3zpg{zpp + Zpg). (15) 

The most constructive way to test the validity of this analysis is to com- 



pare the value of yp at equilibrium, as predicted by iterating eqns.(|T^-|T5|), 



to numerical simulations. In these simulations, the sites are initially ran- 
domly filled with +l's or — I's, so the corresponding initial conditions for 
the iteration procedure are yp\t=o = Zpp\t=o = Zpg\t=o = 0. The results are 
compared in fig.l., where the simulation results compare favourably with the 
approximate analysis, the agreement improving for larger values of p. 

Note that even when p = 1, t/p does not tend to +1, either in the theory 
or in the numerical work. This is because it is possible to have a jammed 
state where, for instance, a g-site surrounded by 4 p-sites may initially start 
at +1 but be unable to change, since if all 4 NN p-sites get fixed into a +1 
state before they have reacted with the central g-site, then the g-site will 
never be able to react and so it will stay as +1 for all time, despite the fact 
that it has Pij = 0. 

IV. INHOMOGENEOUS DESORPTION 

We now turn to an enhanced model studied by Fichthorn, Gulari and 



Ziff ini], which introduces noise into the system in the form of the desorp- 
tion of A and B species from the substrate. They demonstrated numerically, 
later confirmed by mean- field analysis |]T^, that even a small desorption 
rate induces steady-state reactivity onto finite lattices. In our version of the 
model, sites vacated by desorption are refilled by an A or a B as defined by 
the inhomogeneity matrix, which we now call Qij. Qij differs from Pij in 
that now it only applies to sites refilled after desorption - sites vacated after 
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an A + B ^ 2S reaction have an equal chance of being refilled either by an 
A or by a B. Thus, the reaction kinetics alone are the same as the usual ho- 
mogeneous model, and the Uij,Vij and Wij operators without the desorption 
take their simpler form found by setting = ^ in d^-^l). Explicitly, 



8TiUij{S) = A- Sij{Si+ij + Si^ij + Sij+i + Sij-i), (16) 
8t2V,,{S) = l-Si,Si^,„ (17) 
8r2W,,{S) = 1-S,,S,,+,. (18) 

To include inhomogeneous desorption within this formulation, we replace 
Ui, with Ufj, 

Ui = U,, + ^{1 + S,,il - (19) 

where, as in ||^, we introduce a renormalised time scale r and the spin-fiip 
parameter 7, defined by 



1 _ 1 1 1 
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1 - r/rs. 



(20) 
(21) 



The one point spin-correlation rate equation can now be recalculated 
using (I) and ([T6|-[19D, 
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4^^(^^.) = 1^^J{S^J) - 4(1 - ^){{S,,) + (1 - (22) 

This can be solved by using a generating function, G{X,Y,t), defined in 
terms of the time-dependant one-spin correlation function (Sij), 

oo oo 

G(X,r,t)= Y: E X^Y^S,,). (23) 

i=— oo j=—oo 

Combining this with (^) gives rise to a differential equation for G, 

Noting that, except for G(X,Y,t), the right hand side of (^) is inde- 
pendent of time, it is not difficult to derive an explicitly time- dependent 
expression for {Sij) in terms of its initial state, o"jj = (5^^)14=0! 

k,l=—oo 

where Ii{t) is the i*^ order modified Bessel function. In the special case 
^ijQij = it is possible to to rewrite the second term on the right hand side 
of dD as 
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— T — X! ^ '^Qkl) {/j-fc+1 j-l + fi-k-1 j-l + fi-k j-l+l + fi-k j-l-l} ; 

4'''3 k,l=-oo 

(26) 

where for clarity we have introduced 




which obeys the identity 

fi+lj + fi-lj + fij+l + /ij-l 

with the usual Kronecker Delta 
results in an exact expression, 



4 4r 

— fij ^i^djQ 

11 
7 V2r 



7t 
2r 



(28) 



. Substituting ( p8l) into (|25|) and (^) 



(5,,) = 2Q,, - 1 + e-*/^ E (1 - 2Qh + aH)/^-fc (^) • (29) 



kA=—oo 



2tJ ' \2tJ 

So when AijQij = 0, (Sij) — ^ 2(5ij — 1 exponentially as t ^ cx), again in 



agreement with the homogeneous result of |Tl|. With desorption, jamming 
is no longer possible and so now (Sij) 1 when Qij = 1. Although this 
final solution is exact, it is hard to see what physical applications a mixed 
homo/inhomogeneous model such as this one may have. 
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V. INHOMOGENEOUS ONE-SPECIES DESORPTION 



Whilst investigating the monomer-dimer model, Ziff, Gulari and Bar- 
shad [m briefly discussed the additional feature of allowing just the monomers 
to desorb. Physically, this corresponds to the reaction 2C0 + O2 — > 2CO2 
where only the CO can desorb from the substrate, which is a good approxi- 
mation for this reaction at the usual operating temperatures. 

To apply a similar principle to our monomer-monomer model, we extend 
the analysis in section IV to allow for the desorption of A-species only, with 
the inhomogeneity matrix Qij only applying to sites vacated after desorption. 
Thus, the flip-exchange operators are unchanged from (|I^-|T^), but now we 
replace Uij with 



Ui = U., + (1 + S,,) . (30) 

Furthermore, Qij is also taken to be a constant matrix, Qij = q The 



rate equation for the one-spin correlation function (22) is now 



4r|(5,,) = A,,(^,,) - 7(1 - g)(l + (S,,)). (31) 
The definitions of r and 7 have now altered from the previous case. 



1 

T 



1 



1 

T2' 
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(32) 



4r , , 

7 = -. 33) 

T3 



Applying the same generating function ( ]23|) results in a new partial dif- 



ferential equation for G{X, Y, t) acting on an L x L lattice, 



dG G f 1 , 1 , , 1 \ ^} 1-?/^ , t2 



Continuing as before, an explicit time-dependent expression for (Sij) is 
reached. 



L/2 
k,l=~L/2 



2r 



T"3 



t 

I2r 



(35) 



VI. CONCLUSIONS AND DISCUSSION 

We have introduced a methodology for dealing with the effects of gener- 
alised surface disorder on the monomer-monomer reaction process 
A + B ^ 2S, by mapping the system in the reaction-controlled limit onto 
an Ising Model. The two-dimensional rate equations were derived, including 
the very concise one-spin correlation equation (||), and used to study the 
special case of two-site disorder. Here, it was found that the global system 
dynamics are sensitive to the choice of layout of the two different types of 
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site. Catalysts consisting of two different molecules arranged in a regular 
manner, such as on two alternating c(2 x 2) sublattices, allow for no reactive 
equilibrium and will always saturate on finite lattices. Choosing to randomly 
arrange the sites, however, allowing compact clusters of the same site, was 
shown to produce a reactive steady-state. Analysis based on the rate equa- 
tions was used to predict the concentration of A's and B's on the different 
types of site, showing reasonable agreement between theory and simulation 
despite the rather crude approximations involved in the analysis. 

The model was then extended to include desorption from the substrate, 
either by one or both species, and was solved exactly in both cases. 

Extending this work to dimensions other than c? = 2 is straightforward 
once the mapping onto the Ising model has been achieved. Indeed, the rate 
equations for d = 1 can been immediately seen from those given here (p|-[TO|). 
We have focused on d = 2 since the most useful physical application is of 
surface catalysis. 

It should be noted that the definition of inhomogeneity we chose to employ 
here is only one of many ways of modelling surface disorder. For instance, 
requiring that each site be 'hit' a different number of times before adsorbing 
a particle, or assigning a quenched random 'energy' to each site and always 
adsorbing the particles onto the vacant site with the lowest energy, are just 
two alternative possibilities. We intend to study some of these in future work. 
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Figure 1. Plot of p vs. yp\t=oo- The line gives the values predicted by 
the rate equations. Numerical simulation results are plotted as crosses. The 
simulations were performed on a 200 x 200 lattice, and averaged over 100 
runs. 
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